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We present a fast Monte-Carlo algorithm for simulating epitaxial surface growth, based on the 
continuous-time Monte-Carlo algorithm of Bortz, Kalos and Lebowitz. When simulating realis- 
tic growth regimes, much computational time is consumed by the relatively fast dynamics of the 
adatoms. Continuum and continuum-discrete hybrid methods have been developed to approach this 
issue; however in many situations, the density of adatoms is too low to efficiently and accurately 
simulate as a continuum. To solve the problem of fast adatom dynamics, we allow adatoms to take 
larger steps, effectively reducing the number of transitions required. We achieve nearly a factor of 
ten speed up, for growth at moderate temperatures and large D/F. 



PACS numbers: 68.55.Jk, 68. 35. Fx, 81.15.Aa 
INTRODUCTION 

Molecular beam epitaxy (MBE) is a popular technique 
for growing materials, and is also an interesting exam- 
ple of a non-equilibrium statistical process. A simplified 
view of MBE, which captures many of the salient fea- 
tures, is the following. Atoms arrive at a crystal surface 
at a rate F. Once on the surface, the atoms diffuse with 
diffusion constant D until they attach to an existing is- 
land or nucleate a new island. Atoms detach from island 
edges with rate w n — D enp(-ne) where n is the number 
of in-plane nearest neighbors and e is the ratio of effec- 
tive inter-atomic bond energy to the the thermal energy 
(ksT). Over time, adatoms nucleate new islands, islands 
grow by absorbing adatoms, and the surface grows layer 
by layer. 

A common technique for modeling epitaxial growth is 
kinetic Monte-Carlo (KMC). In KMC, atoms are treated 
as individual particles which diffuse across flat substrates 
by taking nearest neighbor hops at a rate AD. Atoms 
which are attached to islands execute nearest neighbor 
hops with a reduced rate (w n ). By representing each 
atom individually, KMC automatically incorporates in- 
ternal noise. In many cases, KMC is a very effective tool 
for simulating epitaxial growth, as it retains its accu- 
racy and flexibility through a wide range of growth envi- 
ronments. However, this dependable accuracy and flexi- 
bility sometimes comes at the expense of computational 
speed. For example, when the adatom density becomes 
large (such as at high temperature) or when the adatom 
diffusion rate is much greater than the island-atom de- 
tachment rate, a KMC simulation spends much of its 
computational effort computing adatom dynamics. To 
compensate for the adatom bottleneck, we devised a new 
algorithm for simulating epitaxial growth, which is faster 
than KMC, while retaining KMC's accuracy and flexi- 
bility. Our method is a KMC approach, but we modify 



the microscopic dynamics of the adatoms to increase ef- 
ficiency, while retaining accurate meso- and macroscopic 
dynamics. The essential step in the technique is to allow 
the adatoms to hop many steps at once. We show that 
this modified dynamics reproduces the true dynamics in 
the regime of interest. 



KINETIC MONTE-CARLO 

The most basic KMC algorithm used to simulate epi- 
taxial growth, which we will refer to as rejection KMC, 
goes as follows. Consider an L x L grid, choose a grid 
point at random, and calculate the hopping rate for the 
atom at that site (w n = exp(— ne)), where n is the num- 
ber of in-plane nearest neighbor bonds. Since adatoms 
have the fastest hopping rates, we let our time-step be 
the time it takes for an adatom to hop. Thus At = 1, 
and the probability for any atom to hop in one time-step 
is w n . We then choose a random number between and 
I, and if that number is less than w n , the atom executes 
a nearest neighbor hop. To complete one time-step, this 
is repeated L 2 times - enough so that on average every 
particle is given one chance to move. 

This algorithm is simple enough to be useful for a wide 
variety of problems in epitaxial growth, including het- 
eroepitaxy. However, it can be quite slow. Since the 
hopping rates differ exponentially, many moves are out- 
right rejected. To avoid this bottleneck, we use the BKL 
algorithm pj, which samples according to the appropri- 
ate rates without rejections. The system we're simulat- 
ing has a very simple rate structure. There are only five 
possible transitions in basic homoepitaxy, corresponding 
to the to 4 possible in-plane bonds. Thus we can use 
the BKL algorithm - often referred to as continuous-time 
KMC. We keep lists of atoms with each bond count - a 
list of adatoms, a list of singly bonded atoms, and list 
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of doubly bonded atoms, and so on. One loop of the 
algorithm involves exactly one transition. The hopping 
rate for each type of atom is w n = cxp(— ne). The total 
transition rate for the system is then W = Ylj—p WjNj 
where Nq is the number of adatoms, N\ the number of 
singly bonded atoms, and so on. The mean waiting time 
for a transition to occur is 1/W. So each loop of the 
algorithm advances time by 1/W (which varies, since Nj 
varies). One of the five types of atoms is then chosen 
with probability 
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and a random member of that set of atoms is hopped. 
In a system with a few fast moving atoms and many 
slow moving ones, BKL runs much faster than rejection 
KMC, and so has become the standard for homoepitaxy 
simulations. It is worth noting that the ideas presented 
in this paper are also applicable to rejection KMC, even 
though we focus on BKL. 



ADATOM DIFFUSION 

As mentioned above, there are five rates in our prob- 
lem. The fastest is adatom hopping (wq — 1), and the 
next fastest is hopping of singly bonded atoms (expo- 
nentially slower at Wi — exp(— e)). So if e = 5 (about 
600-K" for Cu for example), then about 148 adatom dif- 
fusion events occur for every singly bonded edge-atom 
hop - and about 22, 000 adatom events for every doubly- 
bonded edge atom hop. Thus, in regimes where adatom 
diffusion plays a significant role, much of the computa- 
tion is spent calculating the adatom trajectories. 

Our goal is to reduce the disparity in the transition 
rates. We do this by altering what we consider to be an 
adatom event. Normally an adatom transition involves 
one nearest neighbor hop; instead, we allow an adatom 
to diffuse for rid steps. In other words, when we choose 
an adatom, we allow it to execute rid random nearest 
neighbor hops before moving on to the next atom. This 
new n^-step adatom transition must occur at a rate w' = 
1 /rid- The total transition rate is the sum of the rates 
for all the possible transitions 
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which becomes 
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Since the WjS are Boltzmann factors, the total rate is 
dominated by the adatom term (when the temperature 
is low enough). Thus W is often much smaller than W . 



Though the total transition rate has been greatly re- 
duced, the included events have become more complex 
and time-consuming (an n^-step random walk versus a 
one-step nearest neighbor hop). Some increase in effi- 
ciency is realized through this rate reduction technique, 
but there is more to be gained. For the growth conditions 
considered within this paper (high D/F), the step den- 
sity is often quite low, and many adatoms are far from 
an step. Given an adatom on an open terrace (no step 
edges or other adatoms around) the subsequent motion 
is known a priori (in a probabilistic sense). If it is known 
that an adatom is far enough from any obstacles, then 
the adatom can "complete" an n-step random walk in 
one step. The difficulty lies in knowing how much room 
the adatom has. A KMC method that includes long steps 
has been developed in 1-D Q ; however tracking distances 
in 2-dimensions is significantly more challenging. 

What is needed is a way of knowing the distance to 
the nearest attachment site (a site with a free in-plane 
bond). Several methods have been developed to solve 
this issue in various relatedproblems, such as diffusion 
limited aggregation (DLA) [3|. A fast method for doing 
this that has been suggested is to use a hierarchy of 
course-grained grids, where the depth of the hierarchy is 
variable, and determined by the distance from an edge. 
This works nicely in the context of DLA - where the 
growth is irreversible, there is no nucleation, and there 
is only one random walker at a time. But when there 
are multiple random walkers, detachment, edge diffusion 
and nucleation/breakup of dimers, the map hierarchy be- 
comes too cumbersome to be efficient. Methods have also 
been developed to directly compute the distance to an 
edge 0, IE 0- Such methods are very fast when com- 
puting the distance globally; they are too slow, however, 
when only a local update is needed. 

Here, we take a different approach to examining a 
adatom's local environment: we search locally. An n-step 
random walk on square grid potentially covers an area of 
2n(n + 1) + 1 (a diamond of diagonal length 2n + 1). 
In other words, an adatom executing a 25-step random 
walk can end up - or pass through - any of 1301 sites 
(all the sites 25 or fewer hops away). To simulate a 25- 
step random walk, we can compute (a priori) the adatom 
probability density over the 1301-site diamond, and then 
choose the final position of the adatom by sampling ran- 
domly from the probability distribution. This method 
is only accurate if the 1301-site diamond is void of ob- 
structions (steps, attachment sites, and other adatoms). 
Otherwise there would be a sink for the adatom probabil- 
ity density, and the distribution would differ from the a 
prior computed distribution. So to replace a 25-step ran- 
dom walk with a single hop requires searching 1301 sites 
for obstructions, which is computationally cumbersome. 

A simple approximation can reduce the number of sites 
to be searched. Consider that the probability density is 
very low at the perimeter of the 1301-site diamond. In 
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fact, for an n-step random walk, approximately 98% of 
the probability density is contained within a circle of ra- 
dius 2i/n. As a first approximation, we can simulate the 
25-step random walk by searching an area of about 314 
sites, and then sample from the truncated probability dis- 
tribution. That would be faster than searching all 1301 
sites, but still too slow. We could continue to truncate 
the distribution closer and closer to the origin, but the 
approximation becomes less and less accurate. 

Ideally, we would sample from a distribution that is 
as small as possible, yet one that provides for accurate 
dynamics. We start by recognizing that the probability 
of finding a particle at after a 1-step random walk 
is 
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which becomes the diffusion equation 
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in the continuum limit. Given the approximation 
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and retain the same continuum limit. Essentially, we are 
sampling from a distribution which contains only four 
points - (0, ±m), (±m, 0), by allowing an adatom to hop 
a distance m in one of four directions. To make sure that 
the long-hopping adatom doesn't skip over anything, we 
need to search the square that extends out (m-l)-spaces 
- (2m — l) 2 sites. To simulate a 25-step random walk, we 
search a square containing 81 sites, and then hop 5 spaces 
in one of 4 directions. Since each hop of a random walk 
requires the generation of a pseudo-random number, an 
81-site search can be done faster than a 25-step random 
walk. 

Our algorithm is like the standard continuous time al- 
gorithm, but with the following changes. A step distance 
(m) is chosen, so that the associated time step (m 2 ) does 
not exceed the ratio of the adatom hopping rate to the 
next fastest process (usually singly-bonded atoms). All 
the rates w\ — > W4 arc unchanged, but the adatom tran- 
sition rate wq becomes 1/m 2 . Then, if an adatom is 
chosen, we search the box around the adatom to find out 
how far it can hop without hitting something (up to a 
maximum of m). If it has room to hop a distance m, 
then it does, and we move on the the next transition. If 
it has room only for a hop of distance s < m, then it does 
so, but still has m 2 — s 2 time left. So it repeats the pro- 
cess - this time with a maximum step size of V m 2 — s 2 



- until it uses up the full m 2 time-step. That way all the 
adatoms diffuse for the same length of time, but the ones 
near the edges (or other adatoms) take several smaller 
steps instead of one large one. 

Essentially we have replaced the typical adatom tran- 
sition (one nearest neighbor hop) with a rescaled adatom 
transition (m 2 hops), so that the adatoms execute tran- 
sitions with a rate nearer to that of singly bonded atoms. 
Furthermore, we can do the same to the one-bond atoms, 
since they make transitions at a rate much faster than 
the doubly bonded atoms. The general idea is to have 
each type of atom operate on its own time-scale, so that 
the rates become equalized. Of course, for this specific 
problem, there is little reason to go beyond singly bonded 
atoms, so we stop there. It is worth noting that if we were 
to have enhanced edge-diffusion, this multiscale approach 
would provide additional speed-up by also rescaling edge 
diffusion. 



THE ALGORITHM 

We now assemble the ideas presented above into an 
algorithm for growing on a substrate. As in continuous- 
time KMC, each loop of the algorithm is responsible for 
making one transition. The time is then updated accord- 
ing to the total transition rate, and the loop is run until 
we reach the finishing time. 

1. Compute the rates associated with each transition. 
Let the number of atoms with k bonds be N)~; then 
the total transition rate is 
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where m is the hopping distance for an adatom (and 
m 2 the rescaled time-step) and q is the rescaled 
time-step for singly bonded atoms. 

2. Advance time by 1/W. 

3. Choose a type of atom to move. Do this by generat- 
ing a pseudo-random number r between and W. 
If r < F then we add flux (one randomly placed 
atom); if F < r < NqWq do an adatom event; and 
so on. The six possible events are the five atom 
hops (adatom, singly-bonded, etc.) and deposition. 

4. Update the lists of each type of atom, and then 
restart the loop. 

Four of the six possible events are identical to continuous- 
time KMC (adding flux, and hoping atoms with two, 
three or four bonds) . The algorithm for an adatom event 
is as follows: 

1. Choose uniformly from the list of adatoms. 
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2. Scan a square region of side-length 2m — 1. 
Throughout the simulation we keep track of the 
number of bonds available to attach to at a given 
site. Any site with more than zero attachment 
bonds is a potential attachment site. So, we start 
at the innermost square surrounding the adatom, 
and sum up the attachment bonds in all the sites 
in that ring. If that sum is zero, we move on to the 
next concentric square. We keep doing that until 
either the sum is non-zero or we reach the (m— l) th 
ring. Then we know how long our step can be, say 
s. 

3. We then take a step of size s in a randomly chosen 
direction. If s = m, then we are done; otherwise, 
we go back to step 2, but with a maximum step of 
y/m 2 — s 2 . If at some point the adatom attaches, 
we then choose another adatom to finish out the 
time-step. 

The singly-bonded atoms act like the adatoms, except 
that they only take nearest neighbor hops, and only for 
q steps. 

RESULTS 

To arrive at our MSKMC model, we've made some 
approximations regarding the adatom dynamics. While 
the approximations are physically reasonable, it is im- 
portant to verify that the results produced are accurate. 
Since we have only modified the dynamics, any equilib- 
rium results should be unaffected. We can confirm this 
by computing equilibrium island shapes. The shape of 
an island in equilibrium with the surrounding adatoms is 
known exactly by mapping onto the 2-d Ising model 
We took as initial condition a square island on a peri- 
odic grid, and allowed the island to come to equilibrium 
with the surrounding adatoms. The resulting island is 
shown in Fig. ^with the exact (Wulff) solution super- 
imposed. By using a large island (~ 5 x 10 5 atoms, or 
about 250nm for copper), we were able to achieve good 
agreement with the exact solution, without ensemble av- 
eraging. This demonstrates that our method preserves 
the original energetics. 

Since we have modified the dynamics of the adatoms, it 
is important to verify that the new microscopic dynamics 
yield the correct meso- and macroscopic results. A com- 
mon difficulty among continuum and quasi-continuum 
methods is predicting the rate and location for nucle- 
ation of new islands. In models that treat that adatoms 
as a continuum field, the formation of new islands must 
be computed from a mean-field representation of the 
adatom density. In situations where fluctuations play a 
large role, a mean field treatment can give poor results. 
This is especially noticeable at early stages of growth, 
when nucleation - and subsequent breakup - occur very 



rapidly, and correlations in the adatom density can build 
up. MSKMC avoids these problems. To demonstrate 
this, we compare it to ordinary continuous-time Monte- 
Carlo. Both models are allowed to grow starting from a 
flat substrate, and the island sizes are tallied at various 
coverages. The results are shown in Fig. |3 where the 
island size histograms are plotted for coverages of 10%, 
15%, and 20%. There is good agreement between the two 
models. This lends confidence in the accuracy of our new 
adatom dynamics. For this case, the MSKMC simulation 
ran 6.5 times faster than the standard BKL simulation. 

We now test the long-time kinetics, by observing 
mounding. In the presence of an Ehrlich-Schwoebel (ES) 
barrier, a growing surface roughens and forms mounds. 
The ES barrier is a step-edge barrier that reduces the rate 
at which atoms ascend and descend steps. This tends to 
increase the adatom density on the island tops, thus in- 
creasing the nucleation rate on top of existing islands. 
This gives rise to a mounding instability, which has been 
observed in both simulation and experiment 0, 0, 0] . 
The degree of mounding and mound coarsening is sensi- 
tively dependent upon the nucleation on the island tops 
[12j,|ljj, especially when the adatom density is low. This 
is one of the critical challenges of continuum based mod- 
eling; mean-field derived nucleation rates are often inad- 
equate for predicting island-top nucleation 00. The 
sensitivity of mounding on island-top nucleation makes 
this a challenging test for a method which modifies the 
adatom dynamics. 

Rejections are used to implement the ES barrier. The 
ES barrier can be handled using the BKL framework; 
however doing so greatly increases the number of lists re- 
quired, and thus the complexity of the sorting. Since the 
ES barrier is typically small and only arises in certain 
configurations, the number of rejections is low - and so 
rejection can be an efficient alternative to sorting lists. 
All the other barriers in the simulation are handled us- 
ing BKL, and when an ES barrier is encountered, re- 
jection is used. We grew, using both regular BKL and 
MSKMC, 1000 monolayers on a 1024 by 1024 grid. Fig- 
ure |3 shows the surface roughness as a function of time. 
The roughness starts off oscillating; eventually, the os- 
cillations die down, and the overall roughness increases 
due to the mounding instability. The inefficiency of the 
BKL computation makes ensemble averaging over mul- 
tiple realizations prohibitive; however the difference be- 
tween the two results is consistent with observed fluctua- 
tions. Snapshots of the surface at 500 and 1000 monolay- 
ers are shown in Fig. |gj), for both KMC and MSKMC. 
Our MSKMC method ran 7 times faster than the stan- 
dard BKL method. Our method took about 24 hours to 
run the simulation on a 2.8GHz Pentium desktop. 

Mounding can also occur in the absence of an ES bar- 
rier; however it can take a long time to happen. Monte- 
Carlo simulations have shown this effect 15], but only 
when there is enhanced edge diffusion. The increased 
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efficiency of the MSKMC method allows us to simulate 
larger grid-sizes and longer times. In Fig. [5J we show the 
result of growing 1000 monolayers on a 2048 by 2048 grid, 
with no ES barrier and no enhanced edge diffusion. We 
are able to observe decay in the step-density oscillations, 
and also observe the formation of interesting large-scale 
features. This simulation took a little less than two days 
on a 2.8GHz Pentium desktop. 

Epitaxial growth on a vicinal surface is another prob- 
lem of recent interest. In the presence of an ES bar- 
rier, step-flow becomes unstable and the steps can mean- 
der considerably 0. This effect has been demonstrated 
computationally and experimentally [TR Hil Hfl Eol Elj. 
Some systems exhibit meander wavelengths which are 
quite large. For example, Hibino et. al. Hcj show 
experiments involving step-flow on Si, where the mean- 
der wavelength is several microns. Simulating such a 
system requires a very large computational grid. Cou- 
pled with the low temperature and large D/F typical of 
many of these experiments, KMC can be quite slow. Our 
MSKMC method is well suited to computations involving 
large-scale features, low T and high D/F. 

Figure El shows and example of meandering step-flow. 
Both KMC and MSKMC are shown for comparison. The 
computation was done on a 1000 by 1000 grid, with an 
initial step spacing of 20 atoms. The pattern that forms 
is v ery similar to patterns seen for growth on vicinal 
Cu 0. The MSKMC method ran 9 times faster than 
the standard BKL method. For smaller miscut (large 
step-spacing) and lower temperature, the speed advan- 
tage would be even greater. 

DISCUSSION 

There have been several distinct attempts to increase 
the computational efficiency of epitaxy simulation. Here 
we will discuss a couple of these other approaches. 

Absorbing Markov chains 

There are some similarities between our MSKMC 
method and Novotny's MCAMC method [12]. Both 
methods desire to improve computational efficiency by 
rescaling the fast dynamics; however there are signifi- 
cant differences. MCAMC follows the exact dynamics in 
manner that waits for significant changes in the current 
system state before executing a transition. It requires 
knowledge of the transition rates for not only the cur- 
rent state, but also for a number of accessible states. 
There are some situations for which such knowledge is 
efficiently obtained, and as Novotny demonstrates, signif- 
icant speed-up can be seen. Unfortunately, for a general 
epitaxy simulation, it is simply too expensive to compute 
the transition rates for nearby states. It is therefore nec- 



essary to carefully approximate the microscopic kinetics 
in a way that maintains the correct meso- and macro- 
scopic dynamics - which is what we have done. 

Continuum models 

Continuum models, which treat the adatoms as a con- 
tinuous fluid, have been suggested as an alternative to 
KMC |U l^j . Since the adatoms are represented by 
a continuous density, the computational speed is inde- 
pendent of the number of adatoms involved. Continuum 
models with deterministic attachment and detachment 
have had substantial success, but their usefulness is lim- 
ited to growth situations where fluctuations are unim- 
portant. One key source of fluctuations is shot noise. 
Shot noise arises from the fact that atoms attach one 
at a time, as whole particles. Since diffusion-controlled 
growth is inherently unstable, noise can play a vital role 
in determining the shape of islands. An extreme exam- 
ple is the diffusion limited aggregation (DLA) model of 
Witten and Sander 3] . DLA is a diffusive growth model 
with no smoothing processes (such as edge diffusion or 
detachment), which creates branching fractal structures. 
Continuum models lack sufficient noise to generate DLA- 
like clusters. 

There has been some effort to include shot noise fluctu- 
ations in a continuum model [HHl HHH . Most of 
these methods utilize a hybrid approach where adatoms 
are modeled as a continuous density, but islands are com- 
posed of discrete atoms. Once such a ppr oach is quasi- 
continuum Monte-Carlo (QCMC) [H l^. In QCMC, 
atoms are added to island edges as discrete particles, 
with location determined (probabilistically) by the local 
adatom flux into the island edge. This incorporates shot 
noise fluctuations into the continuum formalism. This 
approach has had many successes, including the ability 
to grow DLA-like clusters and reproduce correct thermal 
roughening results |2^, and to some extent model 
multilayer growth with mounding |29j . 

However, these hybrid models still have some short- 
comings. For one, they're not necessarily faster than 
KMC. In a continuum model, the adatom density is de- 
fined over the entire grid, even if there are only a few 
adatoms. If the adatom density is low, then a continuum 
method must take very long time steps when solving the 
diffusion equation in order to be computationally effi- 
cient. This can introduce error into the computation. 
Thus most continuum-based models are inefficient at low 
adatom densities. It is possible for such method to be 
efficient at high densities; however, even that is challeng- 
ing. Since the nucleation (and subsequent breakup) of 
new islands occurs at a higher rate as the adatom den- 
sity increases, the allowable time-step becomes shorter 
as the density goes up. It is a very delicate matter to 
take a sufficiently long time step, while still obtaining 
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acceptable results for nucleation. 

Continuum-based methods bring difficulties beyond ef- 
ficiency. Since the adatoms are replaced by a density, new 
islands must be explicitly created. Often that is done 
using a mean-field estimate of the nucleation rate as a 
function of density. There is a significant literature on 
the subject of nucleation in epitaxy (see reference |30| for 
a review); much of which uses a rate equation approach 
to calculate the dimer formation rate as a function of 
density, temperature, and flux. The results of such calcu- 
lations are used by many continuum-based models. This 
appears to work reasonably well in the absence of any 
step-edge barriers. However, it has been shown yjj 13 
that when there is a step-edge barrier, the rate equation 
approach fails to predict the nucleation rate on island 
tops. Moreover, growth in the presence of a step-edge 
barrier can lead to a mounding instability which is very 
sensitively dependent on the island top nucleation rate. 

Our aim in developing the MSKMC method is to ef- 
ficiently model adatom processes (and other fast pro- 
cesses) without the problems associated with continuum 
methods. Continuum-based methods are very useful 
in situations where fluctuations are less important, and 
when one wants to closely parallel analytical approaches; 
however, when noise and discrete effects become impor- 
tant, a fully discrete method is often better suited. 

CONCLUSIONS 

We have presented an algorithm for simulating epi- 
taxial growth, which is between five and ten times faster 
than continuous-time kinetic Monte Carlo for some inter- 
esting physical regimes. Our multi-scale kinetic Monte- 
Carlo (MSKMC) algorithm is a continuous-time method, 
which rescales the adatom dynamics in both space and 
time. For many realistic growth conditions, the adatoms 
consume most of the computational time. By allowing 
adatoms to take long steps, we reduce the transition rate 
of the dominant process. 

It is worth noting that the relative efficiency of 
MSKMC (as compared to BKL) is dependent on the 
random number generator. Our method replaces ran- 
dom processes (random walks) by deterministic ones (lo- 
cal searching). Due to the computational complexity in- 
volved in generating a pseudo-random number, an effi- 
cient search is faster than a random walk. To be fair, we 
chose to use the fastest reasonable- quality random num- 
ber generator that we were able to obtain (a variation of 
a lagged Fibonacci generator) [il] , which is much faster 
than many compiler default generators. Were we to use a 
slower generator (perhaps to ensure better randomness) 
then our method would be even faster than we've stated, 
relative to BKL. 

To verify the validity of MSKMC, we've computed 
equilibrium island shapes and early-growth island size 



distributions. Both of which match expected results - 
analytical and computational. We were able to repro- 
duce correct equilibrium island shapes without resorting 
to ensemble averaging. Our method also produces step- 
edge barrier induced mounding that is consistent with 
standard KMC. Due to our increased efficiency, we can 
carry out simulations on larger grids and/or for longer 
times, and perhaps observe phenomena previously un- 
seen. 
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FIG. 2: Island size distribution for several coverages, showing 
good agreement between KMC and MSKMC. 




FIG. 3: When a small Ehrlich-Schwoebel barrier is included, 
the oscillations in the roughness diminish over time as the 
total roughness increases. For this case, E/kT — 5, the grid 
is 1024 by 1024, D/F = 10 6 , and the E-S barrier is lkT. The 
difference between the two lines is a result of noise. 




FIG. 4: Surface images at 500, and 1000 monolayers, for 
growth with a small Ehrlich-Schwoebel barrier. 
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FIG. 5: Even in the absence of an ES barrier, the surface 
begins to roughen, and interesting morphology can arise. 




(a) KMC (b) MSKMC 



FIG. 6: Starting from a uniform train of steps, a step-edge 
barrier causes the steps to meander. The step spacing is 20 
atoms, E/kT = 5, E E s/kT = 3, and D/F = 10 6 . 



